#ifndef __CTP23D__
#define __CTP23D__
#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// P2 element, conforming, 3D
// ***********************************************************************
static INT C_T_P2_3D_dof[4] = {1,1,0,0};
// static DOUBLE C_T_P2_3D_nodal_points[30] = 
//               {0.0,0.0,0.0,  1.0,0.0,0.0,  0.0,1.0,0.0, 0.0,0.0,1.0,
//               0.5,0.0,0.0,  0.0,0.5,0.0,  0.0,0.0,0.5,
// 				      0.5,0.5,0.0,  0.5,0.0,0.5,  0.0,0.5,0.5}; 
static INT C_T_P2_3D_Num_Bas = 10;
static INT C_T_P2_3D_Value_Dim = 1;
static INT C_T_P2_3D_Inter_Dim = 1;
static INT C_T_P2_3D_Polydeg = 2;
static bool C_T_P2_3D_IsOnlyDependOnRefCoord = 1;
static INT C_T_P2_3D_Accuracy = 2;
static MAPTYPE C_T_P2_3D_Maptype = Affine;

// base function values
static void C_T_P2_3D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{ 
  values[0] = 1+2*RefCoord[0]*RefCoord[0]+4*RefCoord[0]*RefCoord[1]+4*RefCoord[0]*RefCoord[2]
        -3*RefCoord[0]+2*RefCoord[1]*RefCoord[1]+4*RefCoord[1]*RefCoord[2]
        -3*RefCoord[1]+2*RefCoord[2]*RefCoord[2]-3*RefCoord[2]; 
  values[1] = -RefCoord[0]+2*RefCoord[0]*RefCoord[0];
  values[2] = -RefCoord[1]+2*RefCoord[1]*RefCoord[1];
  values[3] = -RefCoord[2]+2*RefCoord[2]*RefCoord[2];
  values[4] = 4*RefCoord[0]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[0]*RefCoord[0];
  values[5] = 4*RefCoord[1]-4*RefCoord[1]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[1]*RefCoord[1];
  values[6] = 4*RefCoord[2]-4*RefCoord[2]*RefCoord[2]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[1]*RefCoord[2];
  values[7] = 4*RefCoord[0]*RefCoord[1];
  values[8] = 4*RefCoord[0]*RefCoord[2];
  values[9] = 4*RefCoord[1]*RefCoord[2];
}

// base function values
static void C_T_P2_3D_D000(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ 
  // [ 1+2*x*x + 4*x*y + 4*x*z - 3*x + 2*y*y + 4*y*z - 3*y + 2*z*z - 3*z, 
  //   2*x*x - x, 
  //   2*y*y - y, 
  //   2*z*z - z, 
  //   4*x - 4*x*y - 4*x*z - 4*x*x, 
  //   4*y - 4*x*y - 4*y*z - 4*y*y, 
  //   4*z - 4*x*z - 4*y*z - 4*z*z, 
  //   4*x*y, 
  //   4*x*z, 
  //   4*y*z]
  values[0] = 1+2*RefCoord[0]*RefCoord[0]+4*RefCoord[0]*RefCoord[1]+4*RefCoord[0]*RefCoord[2]
        -3*RefCoord[0]+2*RefCoord[1]*RefCoord[1]+4*RefCoord[1]*RefCoord[2]
        -3*RefCoord[1]+2*RefCoord[2]*RefCoord[2]-3*RefCoord[2]; 
  values[1] = -RefCoord[0]+2*RefCoord[0]*RefCoord[0];
  values[2] = -RefCoord[1]+2*RefCoord[1]*RefCoord[1];
  values[3] = -RefCoord[2]+2*RefCoord[2]*RefCoord[2];
  values[4] = 4*RefCoord[0]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[0]*RefCoord[0];
  values[5] = 4*RefCoord[1]-4*RefCoord[1]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[1]*RefCoord[1];
  values[6] = 4*RefCoord[2]-4*RefCoord[2]*RefCoord[2]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[1]*RefCoord[2];
  values[7] = 4*RefCoord[0]*RefCoord[1];
  values[8] = 4*RefCoord[0]*RefCoord[2];
  values[9] = 4*RefCoord[1]*RefCoord[2];
}

// values of the derivatives in RefCoord[0] direction
static void C_T_P2_3D_D100(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ 
  //4*x + 4*y + 4*z - 3, 4*x - 1,       0,       0, 4 - 4*y - 4*z - 8*x,  
  //-4*y,                -4*z, 4*y, 4*z,   0
  values[0] = -3.0+4.0*RefCoord[0]+4.0*RefCoord[1]+4.0*RefCoord[2]; 
  values[1] = -1.0+4.0*RefCoord[0];
  values[2] =  0.0;
  values[3] =  0.0;
  values[4] =  4.0-4.0*RefCoord[2]-4.0*RefCoord[1]-8.0*RefCoord[0];
  values[5] = -4.0*RefCoord[1];
  values[6] = -4.0*RefCoord[2];
  values[7] =  4.0*RefCoord[1];
  values[8] =  4.0*RefCoord[2];
  values[9] =  0.0;
}
// values of the derivatives in RefCoord[1] direction
static void C_T_P2_3D_D010(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ 
  //4*x + 4*y + 4*z - 3,  0, 4*y - 1,       0,                
  //-4*x, 4 - 8*y - 4*z - 4*x,                -4*z, 4*x,   0, 4*z
  //printf("RefCoord: [%f, %f, %f]\n", RefCoord[0], RefCoord[1], RefCoord[2]);
  values[0] = -3.0+4.0*RefCoord[0]+4.0*RefCoord[1]+4.0*RefCoord[2]; 
  values[1] =  0.0;
  values[2] = -1.0+4.0*RefCoord[1];
  values[3] =  0.0;
  values[4] = -4.0*RefCoord[0];
  values[5] =  4.0-4.0*RefCoord[2]-4.0*RefCoord[0]-8.0*RefCoord[1];
  values[6] = -4.0*RefCoord[2];
  values[7] =  4.0*RefCoord[0];
  values[8] =  0.0;
  values[9] =  4.0*RefCoord[2];
}
// values of the derivatives in RefCoord[2] direction
static void C_T_P2_3D_D001(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ 
  //  4*x + 4*y + 4*z - 3,       0,       0, 4*z - 1,    
  //  -4*x,  -4*y, 4 - 4*y - 8*z - 4*x,   0, 4*x, 4*y
  values[0] = -3.0+4.0*RefCoord[0]+4.0*RefCoord[1]+4.0*RefCoord[2]; 
  values[1] =  0.0;
  values[2] =  0.0;
  values[3] = -1.0+4.0*RefCoord[2];
  values[4] = -4.0*RefCoord[0];
  values[5] = -4.0*RefCoord[1];
  values[6] =  4.0-4.0*RefCoord[0]-4.0*RefCoord[1]-8.0*RefCoord[2];
  values[7] =  0.0;
  values[8] =  4.0*RefCoord[0];
  values[9] =  4.0*RefCoord[1];
}
// values of the derivatives in RefCoord[0] direction
static void C_T_P2_3D_D200(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ //[ 4, 4, 0, 0, -8,  0,  0, 0, 0, 0]
  values[0] =  4.0; 
  values[1] =  4.0;
  values[2] =  0.0;
  values[3] =  0.0;
  values[4] = -8.0;
  values[5] =  0.0;
  values[6] =  0.0;
  values[7] =  0.0;
  values[8] =  0.0;
  values[9] =  0.0;
}


// values of the derivatives in RefCoord[0] direction
static void C_T_P2_3D_D020(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ //[ 4, 0, 4, 0,  0, -8,  0, 0, 0, 0]
  values[0] =  4.0; 
  values[1] =  0.0;
  values[2] =  4.0;
  values[3] =  0.0;
  values[4] =  0.0;
  values[5] = -8.0;
  values[6] =  0.0;
  values[7] =  0.0;
  values[8] =  0.0;
  values[9] =  0.0;
}

// values of the derivatives in RefCoord[0] direction
static void C_T_P2_3D_D002(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ //[ 4, 0, 0, 4,  0,  0, -8, 0, 0, 0]
  values[0] =  4.0; 
  values[1] =  0.0;
  values[2] =  0.0;
  values[3] =  4.0;
  values[4] =  0.0;
  values[5] =  0.0;
  values[6] = -8.0;
  values[7] =  0.0;
  values[8] =  0.0;
  values[9] =  0.0;
}

// values of the derivatives in RefCoord[0] direction
static void C_T_P2_3D_D110(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ //[ 4, 0, 0, 0, -4, -4,  0, 4, 0, 0]
  values[0] =  4.0; 
  values[1] =  0.0;
  values[2] =  0.0;
  values[3] =  0.0;
  values[4] = -4.0;
  values[5] = -4.0;
  values[6] =  0.0;
  values[7] =  4.0;
  values[8] =  0.0;
  values[9] =  0.0;
}

// values of the derivatives in RefCoord[0] direction
static void C_T_P2_3D_D101(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ //[ 4, 0, 0, 0, -4,  0, -4, 0, 4, 0]
  values[0] =  4.0; 
  values[1] =  0.0;
  values[2] =  0.0;
  values[3] =  0.0;
  values[4] = -4.0;
  values[5] =  0.0;
  values[6] = -4.0;
  values[7] =  0.0;
  values[8] =  4.0;
  values[9] =  0.0;
}
// values of the derivatives in RefCoord[1] direction
static void C_T_P2_3D_D011(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ //[ 4, 0, 0, 0,  0, -4, -4, 0, 0, 4]
  values[0] =  4.0; 
  values[1] =  0.0;
  values[2] =  0.0;
  values[3] =  0.0;
  values[4] =  0.0;
  values[5] = -4.0;
  values[6] = -4.0;
  values[7] =  0.0;
  values[8] =  0.0;
  values[9] =  4.0;
}
// values of the derivatives in RefCoord[1]-RefCoord[1] direction
static void C_T_P2_3D_Nodal(ELEMENT * elem, FUNCTIONVEC *fun, INT dim, DOUBLE* values)
{
  //printf("Interpolation for P2 element!\n");
    //INT dim =1;
    DOUBLE coord[3];
    INT i;
    //先计算4个节点上的函数值
    for(i=0;i<4;i++)
    {
       // dof on the verts 
        coord[0] = elem->Vert_X[i]; 
        coord[1] = elem->Vert_Y[i];
        coord[2] = elem->Vert_Z[i];
        fun(coord, dim, values+i*dim);
    }
    INT vert0[6] = {0, 0, 0, 1, 1, 2};
    INT vert1[6] = {1, 2, 3, 2, 3, 3};

     for(i=0;i<6;i++)
     {
       coord[0] = 0.5*(elem->Vert_X[vert0[i]]+elem->Vert_X[vert1[i]]); 
       coord[1] = 0.5*(elem->Vert_Y[vert0[i]]+elem->Vert_Y[vert1[i]]); 
       coord[2] = 0.5*(elem->Vert_Z[vert0[i]]+elem->Vert_Z[vert1[i]]); 
       fun(coord, dim, values+(4+i)*dim);
     } 
}
#endif
